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30.0 INTRODUCTION 


The concepts of eigenvalue and eigenvector have played an important role 
in this course, particularly in Unit 13, Systems of Differential Equations, 
Unit 23, The Wave Equation, and Unit 25, Boundary-value Problems. We 
introduced these concepts in Unit 5, Determinants and Eigenvalues. If o is 
a linear transformation from a vector space V to itself, the scalar A is called 
an eigenvalue, and the vector & an eigenvector, if € #0 and of = 2%. In 
this unit we shall assume that V is finite-dimensional and real (typically, R") 
so that this equation can be written as 


Ax = 4x 


where A is a real n x n matrix and x is a one-column n x 1 matrix. We 
shall discuss numerical methods for calculating the eigenvalues 24, 22, . - - 
and the column matrices x,, X2,... representing the corresponding 
eigenvectors, when the square matrix A is given. 


Unlike the problems of our numerical work in earlier units, in which the 
requirements were few and obvious, here we may have various require- 
ments which will influence our choice of method. In some particular con- 
text we may require just the eigenvalue of largest modulus or magnitude; in 
another we may seek the eigenvalue of smallest modulus or the one nearest 
to a particular real number, and in yet another we may want all the eigen- 
values. In each of these cases we may or may not be interested in the cor- 
responding eigenvectors, but in this unit we concentrate on methods that 
give the eigenvector as well as the eigenvalue; both together are called 
an eigensolution. 


When we want just one or two eigensolutions we are likely to use an 
iterative method. If we want the complete eigensystem (all the eigensolu- 
tions) we shall preferably use a method involving similarity transformations 
of A to a simpler similar matrix P~'AP (sub-section 10.1.1 of Unit 10, 
Jordan Normal Form). The iterative methods are treated in Section 30.3 and 
the use of similarity transformations in Section 30.4. 


Following our standard practice in numerical work we must pay close 
attention to possible instabilities: the inherent instability of the problem and 
the induced instabilities of our various methods. We treat inherent instability 
in Section 30.2, and we shall examine each possible induced instability 
whenever we introduce a new method or apply an old method. As so often 
happens in numerical work we shall find that some quite obvious and 
attractive methods can be completely useless for non-trivial problems! 
In particular, the method based on solving the characteristic equation 


det(A — AI) = 0, 


which we found useful in Units 5, 10 and 13 for 2 x 2 and sometimes for 
3 x 3 matrices, is extremely unstable when applied to larger matrices with 
numerical entries. r 


In Unit 8, Numerical Solution of Simultaneous Algebraic Equations, we 
considered possibilities of improving the accuracy of acomputed solution 
by doing just a little more work. In the eigenvalue problem we can do the 
same sort of thing, not only improving an approximate solution but finding 
useful estimates for its error. The relevant analysis is treated in Section 30.5, 
which is optional. 


We shall concentrate exclusively on real symmetric matrices, which occur 
very frequently in practical problems. The mathematical properties of the 
symmetric case are considerably easier than those of the general case, and 
as a corollary of this the computing techniques are easier to understand and 
to analyse. The treatment of unsymmetric matrices is rather more difficult, 
and the only comment we make here is that very few of the results and 
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techniques for the symmetric case carry over unchanged to the general case. 
We shall also assume, though this is not a significant restriction in practical 
problems, that the number of distinct eigenvalues is equal to the order of 
the matrix; that is, all the eigenspaces have dimension 1. 


In Section 30.1 we shall list those mathematical properties of the symmetric 
eigenvalue problem which we shall need for our numerical techniques. 
Some of these properties have been treated in earlier units. 


30.1 MATHEMATICAL PROPERTIES OF THE 
SYMMETRIC EIGENVALUE PROBLEM 


30.1.1 Preliminary Discussion and Results 


(i) What Vector Space to Use? 


In most of this course we have thought of matrices as the numerical 
representation of linear transformations of various finite-dimensional 
vector spaces. In this unit, however, we shall be concerned only with 
matrices, and to avoid having to invent transformations for them to 
represent, it is convenient to work in a vector space V whose elements are 
themselves matrices. We take as our vector space V the set of all one- 
column matrices with n entries, and denote its elements by letters such as 
x, y,.... We shall also need a Euclidean structure on this space: it is given 
by the standard inner product 


x1 "i 
Dp] EN =x to tn 
Xn Dn 


x-y=x'y =y"x (1) 


where x” is the transpose of the matrix x. 


that is 


(ii) The Eigenvector Basis 
Any n x n matrix A defines a linear transformation 
A:x-—— Ax 


in our vector space of n x 1 column matrices. If A is symmetric we can 
apply Theorem 5 of Unit 24, Orthogonal and Symmetric Transformations. 
This tells us that V has an orthonormal basis consisting of eigenvectors of 
A. We shall call the vectors in this basis X4, X2, - - -, X,- The orthonormality 
condition by (1), is 


Xy Xy = XFX, = Örs (2) 


where 6,, means 1 if r = s and Oifr#s. 


(iii) The Modal Matrix 


The matrix whose columns are x,, X2, ---» X, is often called the modal 
matrix, because its columns often represent modes of vibration of some 
physical system (see Unit 13, Systems of Differential Equations, especially 
the television programme). We denote it by 


Because of (2) it satisfies 
X7X=1 


and is therefore an orthogonal matrix (see Theorem 3 of Unit 24). As we 
have verified in sub-section 10.2.1 of Unit 10, X has the property 


AixXay +s An Xia 
AX =) : : =XD 


Axa «++ An Xan 


X7AX=D 8) 


where 


4, 
Equation (3)can also be interpreted by saying that X is a matrix of transition 
that diagonalizes A. 


(iv) Coordinates with respect to an Eigenvector Basis 
Since {x,,..., x,} is a basis, any vector y can be expressed in the form 


Y =X, + °° H On Xy 
with suitable real numbers a,,...,a,, which are the coordinates of y with 
respect to this basis. Since the basis {x,,..., x,} is orthonormal, we have 
yry=yysapt- tay 
As we saw in Unit 16, Euclidean Spaces I (page K278) the coefficients a, can 


be calculated by taking the inner product with x, and using the ortho- 
normality condition (2): this gives 


a, = x7y = y7x, (4) 


(v) Quadratic Forms and Eigenvalues 


In Unit 14, Bilinear and Quadratic Forms, we defined a matrix A to be 
positive definite if and only if the quadratic form y7 Ay is positive for all 
non-zero vectors y. We can express the condition for A to be positive 
definite in terms of its eigenvalues, since 


Y Ay = y7A(@,X, + +++ +0, %,) 
= yT(@AX, + °° + a, A,X) 
=a}, tess +022, (5) 


(from (4)), where x, is the eigenvector of A Corresponding to the eigenvalue 
4,. From Equation (5) we see that A is positive definite if and only if all its 
eigenvalues are positive. 


From Equation (5) we can also derive a lower bound on [Amax| where Armax 
is defined as the eigenvalue of largest magnitude. Equation (5) gives 


ly7Ay| Saill + + alal 
<S (ai + + a2) Amal 
= (YY) masl- 


In particular if y is normalized (has unit length) then we have 


[Y7AY] < |2nax| (6) 
(vi) Measuring Matrices 


Let us denote by M(x), where x is any one-column matrix, the magnitude 
of the entry in x with largest magnitude, so that for example 


leses 


Let us also denote by M(A) the largest of the row sums of the magnitudes 
of the entries of 4: for example, if 


1 2 3 
Az=]2 -1 0 
3 0 - 


the row sums are 
142+3=6, 
24+ |-1|+0=3, 
3+0+ |-2] =5 
and so M(A) = 6. Then it is not hard to show that 
M(Ax) < M(A)M(x). 


This is true, in particular, if x is an eigenvector x,, corresponding to Arex, 
in which case the left-hand side is equal to 


M(emax Xm) = | Amex | Mm) 
We then have that 
|Amax| Mm) < M(A)M (Xn) 
which gives 
|Amex| < M(A) (1) 


By combining this upper bound with the lower bound (6) we obtain the 
useful inequality 


ly7Ay| < M(A), if y7y =1. (8) 


We can also show, by an argument which we do not give in detail here (it 
uses the Schwarz inequality), that this result generalizes to produce the 
inequality 


|y7Ax| < M(A) (9) 


for any pair of vectors y and x for which y’y = x7x = 1, and we shall use 
this result in sub-section 30.2.2. 


In Unit M100 28, Linear Algebra IV, we “ measured” matrices by means of 
column sums rather than row sums. The measure defined here, using row 
sums, is a little more convenient for eigenvalue work. For symmetric 
matrices the two measures are, of course, always equal. 


(vii) Scaling 


There are some theoretical and computational advantages in arranging 
that the entries in the matrix A are of “reasonable” size, neither “too 
large” nor “too small”. This can always be achieved by dividing every 
entry of A by the entry of largest magnitude, so that the magnitude of the 
largest entry in the new matrix is unity. The new matrix has eigenvalues 
which are just those of A “scaled ” by the same factor. We have asked you 
to verify a special case of this in Exercise 4 below. 


Exercises 


1. Find the eigenvalues and eigenvectors, normalized so that xx, = 1, of 
the matrix 


a- 3] 


2. Write down a modal matrix of Exercise 1,and verify that it is orthogonal. 


3. Verify, for Exercise 1, that X~'AX = D in the notation of this sub- 
section. 


4. If Ais the matrix of Exercise 1 write down the scaled matrix for A, as 
described in the text, and calculate the eigenvalues for this scaled 
matrix. Verify that these are the scaled eigenvalues for A. 


LM 30.1.1 


Solutions 


L 


Using techniques from Unit 5 we know that the eigenvalues 
of A are the roots of 


det(A — 2) = 2? -41 -5 =0, 
that is 
à = —1l, à, =5. 


We then solve the systems of equations 


aa] t 


to find corresponding eigenvectors. 


For A, = —1 any eigenvector has the form 


oe] 


where we choose a so that x{x, = 1, that is such that 


1 
@+a?=lora=+—. 


WA 


Then 


x= 


a zbil 
Repeating the calculation for 2, = 5, we find 
OALA fi] ore Ej 
5AL Zl- 
Taking the first of the two alternatives for x, and x, given in 
Solution 1, we obtain for the modal matrix 


“ald 


(There are three other possibilities arising from other choices 
of x, and xz.) We easily check by matrix multiplication that 


eal Jah J-e J 


Since X is orthogonal from Exercise 2, we have x t= XT, 
so that X-'AX = X7AX, and we verify 


areal IB al ali i] 


-1 0 
Bee: 


The scaled matrix is obtained from A by dividing each entry 
by 3, and is therefore 


31 
bay 

The eigenvalues of this matrix are the solutions of 
?-$4-3=0 


which are 4 and — 4. These are the eigenvalues of A divided 
by the same factor 3 that we used to scale A itself. 


30.1.2 The Characteristic Equation and its Instabilities 


In Unit 5 we saw that the eigenvalues of a matrix A are solutions of the 
characteristic equation 


det(A — Al) = 0. 


Having found a solution A, of this characteristic equation, we can then 
solve 


(AA, Dx, =0 


to calculate a corresponding eigenvector. This is the method used in 
Exercise 1 of the previous sub-section. For any non-trivial problem, how- 
ever, the method is quite uneconomic and quite unstable. 


Let us first consider briefly the economy of various methods of finding the 
characteristic equation. Direct expansion of the characteristic polynomial 
det(A — AI) is extremely laborious, as you will easily discover if you try this 
for a matrix of order as small as 6. We might, alternatively, evaluate 
det(A — AJ) for n different values of 2 and solve for the n coefficients 
Qo,..-, @-, the n linear simultaneous equations 


det(A — 4,1) = a, Af + an-1A77' + +++ + a Ay + ao, 
A E T T 


using the fact that a, = (— 1)". One can show that this requires at least 
$n‘ numerical operations. Or we may use the fact that a matrix satisfies its 
own characteristic equation (Cayley-Hamilton Theorem, Unit 5), which 
here means that 


(=1)'A" + a, , A") ++ + aA + ag I= O (the zero matrix). 


We can compute the first n powers of A and equate to zero n entries, say 
the entries in the first row, of the matrix on the left. This method also gives 
a set of linear equations for the coefficients ao,...,@,-,. A matrix 
product involves n° operations, so that the formation of the powers of A 
has something like n* such operations, again a formidable amount of work. 
Next consider the accuracy of this computation. If all the elements of A 
are four-decimal numbers like 0.2679, we are bound to make some com- 
puter-type numerical errors in all the methods we have discussed. In our 
last two methods, moreover, the linear equations turn out to be very ill- 
conditioned, so that even if the order of A is as small as 6, small arithmetic 
errors produce large errors in the solution. In fact there is no known method 
for computing accurately the coefficients of the characteristic equation 
without the time-consuming expedient of using multi-length computer 
arithmetic to reduce the arithmetic errors in the computation. 


Consider next the calculation of the solution of the characteristic equation. 
This is also somewhat laborious, involving iterative methods such as the 
Newton-Raphson method (Unit M100 14, Sequences and Limits II). More 
important, however, is the fact that the calculation of the solution of poly- 
nomial equations is often a very ill-conditioned problem, so that even 
small errors in the coefficients produce large errors in the roots. Consider, 
for example, the polynomial equation of degree 20 whose solutions are the 
integers 1, 2,..., 20 which is 


27° — 210419 +--+ + 20! =0. 


If the coefficient of 4!° is changed by an amount 272°, less than 1077, the 
original solutions 16 and 17 change to the complex pair 


Ayes 217 = 16.780 i281... 


a substantial and most alarming change! 
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The practical use of the characteristic equation therefore exhibits a very 
large amount of induced instability. It is not inherent instability because, 
as we shall see in the next section, the symmetric eigenvalue problem is 
very well-conditioned, small changes in the matrix coefficients giving rise to 
only small changes in the eigenvalues. 


30.1.3 Summary of Section 30.1 


In this section and the Introduction we defined the terms 


eigensystem (page C5) 
modal matrix (page C7) 
scaling of matrices (page C9) 


We introduced the notation 


max (page C8) 
M(x) (page C8) 
M(A) (page C8) 
Main Results of Section 30.1 


(i) The set V of all n x 1 matrices with real entries forms a Euclidean 
space with inner product 


x-y=x’7y, 


Every symmetric n x n matrix A has a set of eigenvectors {x,,..., X,} 
forming an orthonormal basis for V. 


The modal matrix whose columns are x;, ..., X, is orthogonal. 


(ii) The quadratic form yAy is positive definite if and only if all the 
eigenvalues of A are positive. It satisfies 


ly*Ay| < MA) 


if yTy = 1, where m(A) is the largest row sum of the magnitudes of the 
entries of A. In addition, if x7x = 1 


ly7Ax| < M(A). 


(iii) Multiplying all the entries in a matrix by any number (scaling) 
multiplies all its eigenvalues by that number. 


(iv) Unless the matrix A is of order 2 x 2 or perhaps3 x 3,thecalculation 
of the eigenvalues of A by solving the characteristic equation is both 
unstable and uneconomic. 


Technique 


Scale a matrix so that the magnitude of the largest entry of the new matrix 
is 1. 
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30.2 INHERENT INSTABILITIES 
30.2.0 Introduction 


Following our usual plan in the numerical analysis units, we preface our 
discussion of numerical methods by a discussion of inherent instabilities, 
so that we can recognise situations in which accurate solutions are difficult 
or impossible to obtain however good and stable our method might be. We 
therefore examine the effect on the eigensolutions of small errors or 
uncertainties (perturbations) in the entries of A, by looking at the cigen- 
solutions of the matrix A + eB, where £ is a small number and B a 
symmetric matrix. Following the discussion of “scaling” in sub-section 
30.1.1 we can assume, without loss of generality, that all the entries in A 
and B have magnitude less than or equal to 1. 


It is reasonable to suppose that for sufficiently small ¢ the eigenvalue 2,(e) 
of A + eB corresponding to the eigenvalue 4, of A will not differ very much 
from ,. In fact, it can be shown that 2,(e) can be written as a Taylor series 
in powers of £. 


A(e) =4, + k,e + terms in 8, e°,..., (1) 


where the number k, depends on A and B but not on £. We call k,¢ the 
“first-order” perturbation in the eigenvalue 4,, and we hope to find a 
computable expression for k,. We shall not try to estimate the “ higher- 
order” perturbation terms in (1). 


It can also be shown that every entry in x,(e), the eigenvector corresponding 
to A,(e), can be written as a Taylor series in £, having the form 


x,(e) = x, + ez, + terms in e”, e°.., (2) 


where x, is the normalized eigenvector of A corresponding to 4,, and z, isa 
vector depending on A and B but not on £. To fix the length (and sense) of 
the eigenvector we impose the condition 


x,*x,(e) = 1 
so that (2) gives 
X, X, + eX, 2, = 1 6) 


Since this holds for all e, it implies that the first-order perturbation ez, (and 
all higher-order perturbations) in x, is orthogonal to x,. 


From the definition of an eigenvalue problem, we have 


(A + eB)(x, + ez, +*+) = (Ay + ek, + °° Hx + 82, +++) 


Multiplying out the products and writing down only the terms of degree 1 
or less in e, we find 


Ax, + e(Bx, + AZ,) +557 = À, X, + elk, X, + A2) + 77° 


(this uses the rule for multiplying power series, Theorem I-31 on page K665). 

The two power series represent the same function of e and we therefore have 
AX, = 4,%, 

Bx, + Az, = kx, + 4,2, (4) 

The first equation tells us nothing new, but from Equation (4) we can 


proceed to determine the first-order perturbations in the eigenvalue 4, and 
in the corresponding eigenvector x,. 
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30.2.1 The Eigenvalues 


To find k, we have to eliminate the vector z, in Equation (4) of the pre- 
vious sub-section, and this is effected immediately by taking the inner 
product of both sides of this equation with the eigenvector x,. When we do 
this, the term x, - Az, simplifies as follows, since A is symmetric: 


_X,* AZ, = (Ax,) ` Z, 
= Ax, °%, 
= 0; 


and since x, is orthogonal to z,, the Jast term on the right x, - 4,z, is also 0. 
Our formula therefore becomes 


x, Bx, +0 =k,x,-x, +0 
which simplifies to 

k, =x,° Bx, (1) 
since x, is normalized. 
Using Equation (8) of sub-section 30.1.1 we can then write 

lellk,| <]elM(B) <leln (2) 
since the matrix B is scaled and has n entries in each row. 
Now from Equation 1 of sub-section 30.2.0 we have 


12e) — 4,1 < lellk,] (3) 

Hence from Equations (2) and (3) we obtain 
14.) — 2,1 < le | M(B) (4) 
= M(sB) (5) 


Equations (2), (3) and (4) show that small changes, of magnitude at most €, 
in each of the n? entries of A produce “ first-order” changes in the eigen- 
value that are at most only n times as large: the determination of the 
eigenvalues of a symmetric matrix is quite a well-conditioned problem. 


Example | 


k r 1 05 
Consider the matrix A = h. 51 ] 


which we can show has the normalized eigensolutions 
1 fi 

A=15, x =F [i] 

asis aL [ i] 


J2l-1 


and consider the perturbation eB = & fi il 


Our preceding analysis tells us that under such a perturbation the greatest 
possible “ first-order” change in any eigenvalue is ne with n = 2; in fact we 
have here 


ex, Bx, = 7 {l ufi i} Z fi] = 2. 

This is here the exact change; for the characteristic equation of A + eB is 
(l+e-4)?-(0542)?=0 

or (l +e-— 2) = + (0.5 + e) 

giving 2, 42 = (1 + e) + (0.5 + £), 

so that A, = 1.5 + 2e, 2, = 0.5. 


Exercise 


With the A of the previous example, find a symmetric matrix B, all of whose 
entries are either +1 or — 1, for which the eigenvalue of A + £B, cor- 
responding to A; of A, differs from 2, by an amount of modulus approxi- 
mately 2e. Verify that this is the exact change. 


Solution 


The eigenvector corresponding to A, is 


Then 


ex] Bx, = Hell oo +1] [-] 


where, obviously without loss of generality, we have taken the 
off-diagonal terms of B to be +1. If we call the diagonal terms a 
and ĝ, we find 


ex] Bx, = fe(a + $ — 2), 


so that for maximum size we take « = 8 = —1,and the approximate 
change in the eigenvalue has magnitude 2g. 


The characteristic equation of A + eB is then 
(l-—e-A? - (0.5428)? =0 
or 
(l—e-—A4) = (0.5 + 8) 
so that 
Ay, Ag =1—€+(0.5 + 8). 
Then 
A, = 0.5 — 2e, 


and our approximate change is the exact change. 


30.2.2 The Eigenvectors 


To estimate the change in an eigenvector, say x,, produced by the pertur- 
bation £B in the matrix A, we go back to Equation (4) of sub-section 
30.2.0, which is 


Bx, + Az, =k,x, +4,7,. 


This time we want to eliminate k, ; so we take the scalar product with any 
eigenvector x, of A other than x,. This gives 


X,’ Bx, +X, Az, =4,x,°2, (1) 
Since A is symmetric, the second term on the left is equal to 

(Ax,) ` Z, = Às X; ` Z, 
and so Equation (1) can be simplified to give 


x, : Bx, 
‘ned : : (2) 
X,’ 2, Lad 


15 
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This is true for every s # r. The first-order perturbation in the eigenvector 
x, produced by a small perturbation £B in the matrix A is therefore given 
(see Equation 7-46 on page K278) by 


ez = 8 È (x, 2,)X; 


-i Ga) a 
(Gar) 


since x, z, is 0 if s = r and is given by Equation (2) if s # r. 
By Equation (9) of sub-section 30.1.1, the number x,- Bx, has magnitude 
not exceeding M(B), so that the coefficient of x, in (3) has magnitude not 
exceeding 
|e| M(B) = nel 
14, - Al 1A, asl 


since B is scaled. 


(4) 


We see that the coefficient of x, is not large provided that 4, — å, is not 
small. On the other hand, if one or more of the other eigenvalues of A is 
close to 4, then the perturbation in the eigenvector x, is not necessarily 
small in virtue of the small perturbations eB in A. In other words, eigen- 
vectors corresponding to eigenvalues A, and 2, which are nearly equal may be 
ill-conditioned,and the degree of ill-conditioning is proportional to (å, — 2,)™!. 
This result is perhaps not surprising when we look at the limiting case 
when A, and A, become exactly equal. For then there are two orthogonal 
eigenvectors x, and x, with the same eigenvalue A, = A,, and every non- 
Zero vector in the eigenspace <x, , x,> is also an eigenvector. In other words, 
when A, = å, we cannot determine uniquely the directions of the corre- 
sponding eigenvectors, and it is not therefore surprising that we have 
difficulty in determining these directions when A, is close to A,. 


Example 
The matrix 
1 0.01 
4 = [aor 1 | 


has the normalized eigensolutions 
= 1 1) _ fo7 
a= toy» Je] = [83], 
= —_17 1 0.71 
I= 099,41 = [_ Hl z Bes 


We consider the effect of a perturbation eB, with e = 0.04 and 


1 0.5] - 
a=) a 


and to this end we calculate the eigensolutions of the matrix 


L04 0.03 
AERAR Fe nee 


The eigenvalues are 1.05 and 0.95, not very different from those of A in 
accordance with our expectations. 


The normalized eigenvectors are multiples of the vectors 


LM 30.2.2 


and they are clearly very different from those of A. In fact if we take the 


3 3 
first eigenvector to be GHE and choose k, as suggested in sub-section 


30.2.0, so that we can write 


Galla 


with the coefficient of x, equal to 1, we find 


and 


This means that the perturbed eigenvector corresponding to x, has ac- 
quired a large component, of magnitude 4, in the direction of x! 


You should realize that the formulas (1) of sub-section 30.2.1 and (3) of 
this sub-section do not here give very accurate expressions for the pertur- 
bations. This is because the number 2, — A, = 0.02 is small even compared 
with our particular e, so that “higher-order” effects are here important. 
The general conclusions, of course, are not affected. Even “ first-order” 
perturbations are large in eigenvectors corresponding to close eigenvalues, 
and this is enough to guarantee the ill-conditioning of these eigenvectors. 
Higher-order effects, on the other hand, cannot affect our conclusion about 
the good conditioning of the eigenvalues. In fact we can show, by more 
advanced techniques which we shall not give here, that Equation (2) of 
sub-section 30.2.1, giving a bound for the first-order perturbation of an 
eigenvalue, is true without the qualification “first-order”. In other words, 
for any perturbation B, however large, no eigenvalue of A + £B can differ 
from some eigenvalue of A by more than M(B) in magnitude! 


Exercises 


1. In the example of this sub-section, show that the perturbation in the 
eigenvalue A,, obtained from Formula (1) of sub-section 30.2.1 is 
smaller than the true perturbation, but that the latter still satisfies 
Equation (2) of sub-section 30.2.1 and verifies the last remark of the 
text. 


2. If the perturbation of the example of this sub-section were applied to 
the matrix of Exercise 1 of sub-section 30.1.1, why is the effect very 
small on both eigensolutions ? 


Solutions } 
1. The first-order perturbation in å, is 


ex1 Bx, = 0047 u nfo | Z fi] =0.02, 


whereas the exact change in A, is 0.04. This exact change is 
smaller in magnitude than eM(B), which is 0.04(1.5) = 0.06. 


2. The eigenvalues of the given matrix are 2; = —1, 4, =5, 
which are well separated. Hence the eigenvectors, as well as 
the eigenvalues, are well-conditioned, affected only slightly by 
a small perturbation in the matrix. 
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30.2.3 Summary of Section 30.2 


In this section we defined the terms 


perturbation of a matrix l (page C13) 
first-order perturbation of an eigenvalue (page C13) 
first-order perturbation of an eigenvector (page C13) 


We introduced the notation 
Afs) (page C13) 
x,(e) (page C13) 
Main Results of Section 30.2 


(i) The first-order perturbation ek, of the eigenvalue 2, of a symmetric 
n x n matrix A caused by the small perturbation €B of A is bounded as 
follows: 

lellk,| < le] M(B) < lel” 


In other words, the determination of the eigenvalues of a symmetric 
matrix is quite a well-conditioned problem. 
(ii) The first-order perturbation ez, of the eigenvector x, of A satisfies 


eed Fee i) * 


and 
lellx, : Bx,1 _ lel M(B) nle] 
-Aal SA-A] 


In other words, eigenvectors corresponding to eigenvalues A, and 4, which 
are nearly equal may be ill-conditioned. 


ee 


30.3 ITERATIVE METHODS FOR INDIVIDUAL 
EIGENSOLUTIONS 


30.3.0 Introduction 


This section and Section 30.4 present useful practical methods for computing 
eigensolutions. In this section we examine iterative methods, which are 
most convenient if we want just one, or perhaps a very few, of the eigen- 
solutions. We shall discuss two important iterative methods. The first is 
called direct iteration and this, as we shall see, can produce only two 
eigensolutions, those for which the corresponding eigenvalues are farthest 
from zero in the positive or negative directions. The second is called 
inverse iteration, and this can produce any eigensolution and is therefore, in 
some senses, a more powerful method. 


30.3.1 Direct Iteration: Theory 


We have already observed (in sub-section 30.1.1) that any vector y in V 
can be expressed as a linear combination of eigenvectors in the form 


Yy =X, HAX +°°' +4,X, (1) 


What we try to do with iterative methods is to operate on y in such a way 
that one of the coefficients a, is magnified at the expense of all the others, so 
that we transform y into a good approximation to some eigenvector. 


The simplest operation uses the linear transformation A itself. Weconsidered 
the effect of applying such a transformation repeatedly in Unit 5, sub- 
section 5.3.2, but a little generalization produces more useful techniques. 
We therefore consider a linear transformation with the matrix A — pl, 
where p is some real number, and observe that 


(A — pl)x, = AX; — pX; = (Ar — BX» (2) 


where x, and A, represent an eigensolution of A. It follows that the cor- 
responding eigensolution of A —pl is 


x, and 2, — p. 
From Equations (1) and (2) we find 

(A — pl)y = ay (Ay — p)Xy + 0222 — p)X2 + 71+ + An — PIX 

(A = pl)?y = a, (å, — p)?xy + alz — p)?X2 +77 + aln — P)’ Xn 
and in general 


(A — ply = a, (Ay — P)'X, + a(z — p)X2 + ` + (Ay — ie 
Now as r increases the right-hand side of (3) is increasingly “ dominated ” 


by the term a,(A, — p)"x,, for which | A, — p| is the largest of all the |4; — p|, 
(i = 1,..., n) provided only that a, # Oand that the k so defined is unique. 


We can then write (3) in the form 


A, -PY ; 
(4 plyy = Oa — PY {a (z EV a boo bo 
xP. 


and as r increases every term in the bracket other than a, x, gets smaller 
and smaller, and (A — pl)'y becomes more and more nearly proportional 
to the eigenvector x,. (We saw an example of this in the last sub-section 
of Unit 5.) 
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The same calculation also gives us the eigenvalue A, . For sufficiently large 
r we have 


(A — ply*'y = (ay — py * {as % + small vectors} 
(A — pi'y = (Ay — p)'{ax X + small vectors}, 


and it follows that as r gets larger and larger the ratio of any corresponding 
non-zero entries in the vectors (A — pl)’*'y and (A — pl)'y gets closer and 
closer to A, — p. With a prescribed p we therefore compute a close approxi- 
mation to the eigenvalue ,. 


This method of computing an eigensolution is called direct iteration with 
A — pl. Two important questions are: how do we choose the number p, 
and what particular eigenvalues A, can we actually compute by this method ? 
To answer these questions consider the eigenvalues located at points on 
the real line, shown below. 


o i H 
4+ isP |e] a PI} 


Then the values of |A, — p| are the distances of these points from the point 
pon the real line. If p is closer to the largest eigenvalue 2, than to the small- 
est Ay, then the largest value of |à, — p| is |As — p|. If p is closer to Ag, 
then the largest |2, — p| is |24 — pl. It is clear that the only eigensolutions 
to which direct iteration will converge are those with eigenvalues 2, and 
Ag, that is to the eigensolutions corresponding to the two extreme eigen- 
values. 


This answers both our questions, though of course a decision about what p 
to choose in order to guarantee convergence to A, say, rather than to Ag, 
will need either some prior knowledge of the approximate values of 2, or 
Ag, or a process of trial and error. 


But there is another interesting thing about the choice of p which concerns 
the rate of convergence of our iterative process. To examine this let us 
change our notation, labelling the eigenvalues so that 4, > 2, > ++ > Ap, 
as shown below. 


—- i sg 
An Anat An M 


Let us suppose that we want to calculate 4,, the smallest eigenvalue. The 
relevant equation, from (4), is 


(A ~ PI = On -pY fax, +a, Aa Jsu 


For convergence, as we have seen, we must choose p so that 
lå =pl < [4 — pl; 


this ensures that |2, — p| is the largest of the |A, — P|. Provided «, ¥ 0. 
the rate of convergence clearly depends on the rates at which the various 
|4:—pl/|4, — p| get smaller as r increases, and we would like to choose p 
so that the /argest of these numbers is as small as possible. 
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Now the “ worst offenders” are clearly 2, and 2,_,, and the best choice for 
p is that for which these give equal and opposite ratios, that is 


-P An-i =P 
Laps hap ot A,—p=—-(,-1—p), 
giving 
P = $à, +4,-1). 
Any other choice would make one of these ratios larger than it need be! 


By precisely the same argument we can show that the best choice of p for 
convergence to A,, the largest eigenvalue, is 


P =}, + 23). 


Example 

Suppose that we have a matrix of order 3, with eigenvalues 5, 3 and — 1. 
Taking p = 0, that is iterating directly with A, we converge to the eigen- 
solution with À = 5, and the ratios which govern the rate of convergence are 
ł and 3. The choice 


p=33-N=1 
gives 
A, —p=4,4, -P =2, 1, -p= -2, 


and we still converge to the largest eigenvalue but at a rate governed by the 
ratios 


22-P 
2 =P 


23 -P 
A,-p , 

both of which are 4. Since (4)’ decreases considerably faster than (ły, as 
r increases, we have achieved a much better rate of convergence. 


It is only fair to remark here that a good choice of p requires even more 
knowledge about the distribution of the eigenvalues. But trial and error can 
still play its part, and in early work with digital computers this was quite a 
powerful method, since p could be varied at will and the effects of conver- 
gence could be observed on something like a television screen. 


Exercises 


1. A certain matrix has eigenvalues — 5, 0, 2 and 5. For what values of p 
would we, by direct iteration with A — pI, obtain fastest convergence 
(i) to the largest eigenvalue, (ii) to the smallest? What, in each case, are 
the ratios governing the rate of convergence? 


2. Can you see how you might get both the largest and smallest eigen- 
values, and the corresponding eigenvectors, almost simultaneously by 
direct iteration with A(that is, p = 0), in this particular case when 
aA,=—-A,? 


Solutions 


1. Let us represent the eigenvalues diagrammatically. 


Me Ag Az M 
_— ———> ———o— e 
-5 o 2 5 


(i) For convergence to the largest eigenvalue, which is A, = 5, 
we choose p to be the average of the largest and smallest 
unwanted eigenvalues, which is 


H-54+2)=-2. 
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i 
The ratios ri E 


governing the rate of convergence are 


1 
7. ay ry representing a reasonable rate of convergence, 


(ii) For convergence to the smallest eigenvalue, which is 
A, = —5, we choose p to be the average of 0 and 5, which 
is į. The relevant ratios are 1%, Tz, 73, and we get faster 
convergence than in the previous case. 


2. Direct iteration with A will give, for sufficiently large r, 
Ay = U{(a,x, + (—1)'a, x,) + negligible vectors} 


So 42, giving A, and —A, = ,, is the ratio of corresponding 
entries of A’*”y and A’y. Having obtained 2,, we can then 
write i 


P, 


=a X; +a,x, 


= OX, F nX 


from which we easily find multiples of the required vectors x, 
and x,. 


30.3.2 Direct Iteration: Practice 


It is quite uneconomic to compute successive powers of the matrix A — p1, 
a procedure which appeared to be needed in our discussion of the theory of 
direct iteration. We prefer to express the method in terms of an iterative 
algorithm of the form 


yo") = (A — ply, y arbitrary, 

which you can easily see is identical, at least in theory, with the statement 
yor) =(A —ply* 1y(0) 

implied in the previous sub-section. 


In practice we go just a little further, in the attempt, already mentioned in 
sub-section 30.1.1 (vii), to keep all our numbers of reasonable size. We do 
this by dividing each successive vector by its entry of largest magnitude, so 
that the new vector has the magnitude of its largest entry equal to unity. 
This produces the algorithm for direct iteration in the form 


y arbitrary (with largest entry = 1) 
zt) a5 (A — phy” 
yer) aa z+ Dmt +D), (1) 


where m(z) means the entry in z with largest magnitude. (Note that 
Im(z)| = M(2).) 


Apart from keeping the numbers to reasonable size this device has two 
other advantages. First, »(z*') is the best current estimate of the number 
2, — por A, — p, to which we are converging, at least for an r large enough 
to ensure that successive vectors have the largest numbers in the same 
entry. Second, and more important, this process answers another question : 
what happens if the initial guess y is orthogonal to the eigenvector we are 
hoping to find, that is, if «, (or a) = 0 in the equation 


à 
y= Jax, 2 
r=1 


In theory, and certainly with exact arithmetic, we could converge to some 
other solution! In practice we want to be sure that we have converged to one 
of the extreme solutions, and our scaling process virtually ensures this 
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For with computer arithmetic the division z“ * fm(z"* 1) will sooner or 
later (and usually sooner!) involve rounding errors, so that the computed 
y*” will equal the true y"*” plus a vector of rounding errors. This will 
include a small multiple of the required x,(or x,), and this multiple~is 
steadily magnified as r increases. Convergence will be slow at first but 
ultimately certain, and here rounding errors act to our advantage! 


Example 


2 1 3 
For the matrix A =|1 —1 1 
3 1 4 


we Start with 
yOu} 1 17, 


and (taking p =0) produce the following results in three-digit floating 
point arithmetic. The column labelled 4, is m(z‘) and is the current esti- 
mate of the dominant eigenvalue. (You need check only the first one or two 
steps.) 


r zT yor LA 
ot 1 ) u 1 1] 

1 1 8) [0.750 0.125 1} 8 
2 [4.63 1.63 6.38) [0.726 0.255 1] 6.38 
3 (4.71 1.47 6.43] [0.733 0.229 1) 6.43 
4 [4.70 1.50 6.43] [0.731 0.233 1] 6.43 
5 [4.70 1.50 6.43) [0.731 0.233 1] 6.43 


To three-figure accuracy we have 
2°) = Ay 

and 
y = 259/643 = y. 


To this accuracy 6.43 and y‘®’ is an eigensolution, and we can confidently 
expect that it is the dominant eigensolution. 


Exercises 


1. (i) Given that the eigenvalues of A in the last example are approxi- 
mately 6.4, — 1.3 and —0.1, what value of p would you use for 
fastest convergence to the dominant eigensolution by iterating 
with A — pl. 

(ii) Taking p = —1 for simplicity, perform two steps of the process, 
using three-digit arithmetic and starting as before with y% = 
0 1 1F. 

(iii) What is the best estimate of the dominant eigenvalue of A given 
by this calculation? 


2. The matrix 1 
13 20 04 
A =|2.0 19 -16 
04 -16 3.1 
has eigenvalues 2, = 4.5, A, = 2.7, å, = —0.9. An eigenvector corres- 
ponding to A, is 
x, =(0.5 1.0 —1.07. 
(i) Show that the vector 
y=(04 08 1.0]7 


(ii) 


is orthogonal to the eigenvector x,, so that a, = Oin 


3 
= Fax. 
r= 
Show, however, that if we usc three-decimal floating-point 
arithmetic the iterative method of this sub-section, with p =0, 
will converge to the dominant eigensolution 2, and x,. 


Solutions 
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1. (i) The best value of p is half-way between — 1.3 and —0.1, 
so that p = —0.7. r 
(ii) Taking p = —1, we have 
3 1 3 
A-pl=|1 0 1 
Je) (a) 


and two steps of the iteration give 


r zT yor A+ 
of 1 ou 1 1] 
1 Toa 9) [0.778 0.222 1) 9 


2 [5.55 1.78 7.55] [0.735 0.236 1} 7.55 


(iii) Atthis stage the vector appears to be converging nicely, 
and the current best estimate for the required eigenvalue 
isà — p =7.55. 


With p = —1 this gives A = 6.55. We are converging, as we 
would expect, just a little faster than we did in the previous 
example. 


2. (i) WithxT=[0.5 1.0 —1.0],y® = [0.4 0.8 1.0]7,the 


inner product 
x] y = (0.5)(0.4) + (1.0)(0.8) + (— 1.0)(1.0) =0 


Hence y is orthogonal to x,, and «, = 0 in 


3 
y= $ ax. 


r= 


(ii) With y as first approximation the iteration proceeds 
as follows: 


% 1.3 20 04 
A =|2.0 19 —1.6 
0.4 —1.6 3.1 


r 20T yor 


O (04 08 10) [04 08 10 ] 
1 [2.52 0.72 1.98] [1 0.286 0.786] 


There are already rounding errors in the formation of yo) 
from z", so that a multiple of x, is now present in y”, and 
we shall converge to x, since 2, = 4.5 is the largest eigenvalue 
and we are taking p = 0. (On the O.U. computer we get the 
answer correct to three decimals after 56 iterations.) 


30.3.3 Direct Iteration: Induced Instabilities 


Direct iteration for the extreme eigensolutions has virtually no induced 
instability. It will always produce reasonably accurate eigenvalues, and the 
accuracy of the computed eigenvectors will depend, as one would expect, 
on the inherent instability of the problem. We can see this with the help of 
a backward error analysis, similar in spirit to the one we used in Unit 8 
for the Gauss elimination method. 
The iteration method with p = 0 is described theoretically by the equations 

2041) = gy 

yor) = 20* Dm(g +9) (1) 
that is, by 

Ket Dyle+ 1) = Ay 
where 

ket = m(t +) 


Owing to rounding errors in the computation of Ay, the calculation will 
actually give 


ert Dylr #1) = Ay 4. go (2) 
where f is a vector of rounding errors. 
Now we would like to be able to write this equation in the form 

Kort dyer +1) =(A+ ZMH (3) 


so that we could say that the vectors we have actually produced at this 
stage were obtained from the perturbed matrix A + Z” rather than from 
A itself. 


We would then like some estimate of this perturbing matrix Z“. Now the 
vector y“ has been normalized so that its largest entry is 1. Suppose that 
this largest entry is the sth, and let 

e=[0 0-0 1 0---0)7 


be the vector with zeros everywhere except for a 1 in the sth position, so 
that e, is in fact the sth column of the unit matrix J. Now write 


0---0 FAK 0---0 
Z =e] = | Sa i | 
0-0 AG 0-0 
This is the matrix whose columns are all zero except the sth, which is £ 
Then, since the sth entry in the vector y is 1, we have 
Zoy = 1, 


and from Equations (2) and (3) we see that the matrix Z" of (4) satisfies our 
requirement. 


(4) 


We can now estimate the accuracy of the eigensolution computed by the 
iterative method, in which we terminate the process when two successive 
vectors y“? and y”*” are the same to the number of figures to which we 
we are working. By (3), the resulting vector y**” and the number k”*”? 
constitute an eigensolution of A + Z“?, and we can then apply the theory 
of Section 30.2. For the eigenvalue, Equation (5) of sub-section 30.2.1 gives 
the result 
|ko*)) — trueeigenvalue] < M(Z™) 
= Mf), 


from Equation (4), and this is the magnitude of the entry of f°” of largest 
magnitude. 
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Now if no entry of A exceeds unity in absolute value, and since no entry of 
y exceeds unity in absolute value because of (1), we see that the entries of 
f can hardly exceed 


me 
nx27', 


where £ is the number of binary digits to which our computer performs its 
arithmetic. Indeed this is guaranteed if we exercise some slight extra care 
with the arithmetic, as described on page 35 of Unit 8. 


The computed vector, of course, could be far less accurate, with a possible 
error of 


nx27-t 


ay, 


in its components, and this could be large if any A, is near to 2,. But this 
merely reflects the possible inherent instability of the eigenvectors, and we 
can hardly blame our method for this! 


Induced instability will occur, however, if we try to find any other solution 
by direct iteration starting with a vector y which is orthogonal to the 
dominant eigenvector, that is with 


YO = ag Xp + aX H + OX, 
{assuming that A, is the eigenvalue of largest magnitude). 


If we have computed 4, and x, exactly, then we can find such a y® by 
the Gram-Schmidt method described in Unit 16, Euclidean Spaces I, 
Starting with some arbitrary vector y, and taking 


yO =y- ax, 
with 
a = xy if x7x, =1. 


In theory direct iteration with A will now Converge to a different eigen- 
solution; but in practice the method will fail because (a) we will not have 
X, exactly, (b) we cannot compute æ, or y exactly, and (c) even if we 
could, the rounding errors in the iterative Process, which worked to our 
advantage in computing the dominant eigensolutions, now magnify as 
the calculation proceeds, and insist on Producing once again one of 
these dominant solutions. 

Exercise 


For the matrix lo i 2 '] the dominant eigensolution is i 


à =1.01,x,= [1 _1]7 
scaled so that M(x,) = 1. 
Using direct iteration with the starting vector 
yO = [0.7 17 
and using one-digit arithmetic, show that we obtain an accurate eigen- 
value but an inaccurate eigenvector. Explain this result. 
Solution 


The iteration produces 


r yor yor à 


o [0.7 1) [0.7 1) 
1 [0.7 1.0) [0.7 1.0) 10 
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The method converges rapidly (in one step!) and the computed 
eigenvalue is very accurate, but the computed vector is rather 
inaccurate. This occurs because the eigenvalues are 1.01 and 0.99, 
very close together, and the rounding errors of the arithmetic are 
equivalent to the solution of a slightly perturbed matrix for which 


(both) eigenvectors may be quite different from those of the 
original matrix. 


30.3.4 Inverse Iteration 


We now discuss a method of iteration which will give any required eigen- 
solution. A hint of the idea is given by the title inverse iteration. In fact, 
inverse iteration with the matrix A — pI, where p is any number, is effect- 
ively direct iteration with its inverse (A — pl)~?. 


For the theory of inverse iteration we need the eigenvectors and eigen- 
values of A — pI in terms of those of A. We start from the formula 


(A — pl)x, = (2, — p)x, 
which we have already used for direct iteration with A — pI. Provided 
4, # p, the matrix A — pI is non-singular, and we can therefore pre- 
multiply both sides by (2, — p)“ + (A — pl)7} and obtain 

(A, — p)7*x, = (A — pl)™'x,. 
That is, if 2, # p, then x, is an eigenvector of (A — pl)~' corresponding 
to the eigenvalue (4, — p)~!. 
The effect of operating on any vector 

y =a,xX, +°°° +4,X, 
repeatedly with (A — pl)~' is then given by the formula 

(A — pl) fy = @y(A, — p) X +: + aln — P) Xa 


which has obvious analogies with Formula (3) of sub-section 30.3.1. As r 
increases, the term in the sum that increases fastest is the one for which 
1, —p)7!| is the largest. This, of course, is the term for which |4, — p| 
is the smallest. As r increases, therefore, the right-hand side gets closer 
and closer to a multiple of the eigenvector x, whose eigenvalue is closest 
to p. 

The corresponding iterative algorithm, again, can be written in a form 
analogous to that of Equations (1) of sub-section 30.3.2 relevant to direct 
iteration, and would satisfy the equations: 


y arbitrary (with largest entry = 1) 
2041) = (A — ply ty (1) 
yt 1) = z” + Dim(2e* D) 
N 
Then y” gets nearer to an eigenvector x, as r increases, the corresponding 
eigenvalue A, is the one nearest to p, and its computation comes from the 
equation 
mz” 40) = (A — p's 
the approximation getting steadily more accurate with increasing r. 


The choice p = 0, for example, produces the eigenvalue of smallest mag- 
nitude, and this is often required in practical problems. We can also solve 
immediately the problem of finding the eigenvalue nearest to a particular 
given value of p. It is in fact clear that inverse iteration, which can find 
any solution with suitable choice of p, is far more powerful than direct 
iteration which, whatever choice we make for p, can only produce the 
algebraically largest and smallest eigenvalues. 
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In practice we change the algorithm (1) in one important way. It would 
appear that we have to invert the matrix (A — pl). and we have seen in 
Unit 8 that this is a somewhat laborious operation. It is quicker to solve 
linear equations, and we therefore write the second of Equations (1) in 
the form 


(A = pN +? ay 
which is theoretically the same but computationally quite different! 


Notice, in particular, that for each r the linear equations have the same 
matrix A — pl but different right-hand sides. In Unit 8 we saw that the 
largest part of the Gaussian elimination process with interchanges is the 
reduction of the relevant row-permutation of A — pI to an upper triangle 
U, and this need be done once only. To calculate each new z“*™ we use 
the same multipliers and interchanges on y“ as we used in reducing 
A — pl to U. This converts y to c® such that 


U2lt +) = ol), 


a system of equations from which z*") is found by back substitution. 
This is far more economical than computing the inverse (A — pl)~' and 
then evaluating (A — pI)“ ty” to get 2¢*!). 


Example 


Let us iterate inversely with the matrix 


The first step is to perform the relevant triangular decomposition, and 
since in Unit 8 we only sketched the decomposition process corresponding 
to interchanges in Gaussian elimination we will go through it here in some 
detail. We proceed to perform the elimination with interchanges, but 
keeping close account of the positions of the original rows in A. The 
numbers in brackets refer to the rows of A, and we start with 


7 ES 4 (1) 
1 -1 -4 (2) 
4 -4 -8 (3). 
We interchange rows (1) and (3) and then perform the elimination. The 


multipliers are written on the left, and for simplicity we use exact 
arithmetic. 


4 -4 -8 3) 
-4 1-1 -4 (2) 
-4 2 1 4 (1) 
4 —4 -8 (3) 
0.0 -2 (2) 
0 °3 8 (1) 


We now interchange the second and third rows of this matrix, and then 
eliminate, obtaining 


4-4 -8 O 


0 0 oO -2 (2). 


This shows that the upper triangle U comes from A, which is A with rows 
written in order 3,1, 2. We know that the non-diagonal elements of the 
L matrix (in A = LU) are the negatives of the multipliers, but what are 
their positions in the L matrix? Well. we just look for the multipliers 
associated with the rows of 4. The first is the third row of A. and this has 


no multipliers. The second is the first row of A. and this has the single 
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multiplier =} The third is the second row of A, and this has multipliers 
—t and 0 in this order. The required L matrix is therefore 


100 
| | 0 
401 


and we easily verify that 


A L U 


4 -4 -8 1 0 f4 -4 -8 
2 1 4/=]4 1 offo 3 8 
1-1 -4] [4 0 ajlo o - 


Now we start the inverse iteration, taking the arbitrary y™ as 
yOr = [0.32 1.0 —0.32]. 


The co corresponding to y is obtained by subjecting y to the same 
row interchanges and multipliers that took A to A and thence to U. In 
other words, 


Le = [-0.32 0.32 1.0)’, 

and solving these linear equations by forward substitution we find 
c = (—0.32 0.48 1.08). 

Our z” is then obtained by back substitution in the equations 
U2) = c0, 

and we easily find 
z§) = —0.54, z$} = 1.60, z{) = 0.44. 

Then 
yT = [0.275 1 —0.3375] 

and the current approximation to the eigenvalue is obtained from 
m(z)=147} (since p = 0), 

giving 
å = 1/1.60 = 0.625. 


Compared with direct iteration, the only extra feature of the inverse 
iteration method which could possibly give rise to induced instability is 
the solution of the relevant linear equations. This we know to be a stable 
process if we use Gauss elimination with interchanges. But however 
stable the method, we may have an ill-conditioned mathematica! problem 
if the selected p, in inverse iteration with A — pl, is very nearly an eigen- 
value of A. For then the matrix A — pl is almost singular, and as we saw 
in Unit 8 this is just the situation in which ill-conditioning is most manifest. 


It is extremely unlikely, of course, that our selected p will be very near to 
an eigenvalue, and even then, as we shall show in sub-section 30.4.3, there 
is no trouble if the vector y’, say, in Equation (1), is not nearly orthogonal 
to the eigenvector x, corresponding to the eigenvalue 2, which is near to p. 
As we saw with direct iteration, for any significant value of r, that is after 
a fair number of steps, the vector y” will have picked up a good compo- 
nent of x,, if only through the magnification of earlier rounding errors, 
and we have no induced instability. 


Exercise 


2 13 
The eigenvalues of the matrix A = f -1 j in the example of sub- 
3 14 


section 30.3.2 are approximately 6.4, — 1.3 and —0.1. If we iterate inversely 
with A — pl, with p = 4, to what eigensolution shall we converge, and 
what are the main factors governing the rate of convergence? 

Solution 


The number p = 4 is closest to the eigenvalue 6.4, so that we 
shall converge to the dominant eigensolution. We have 


6.4 — p =24, -1.3 — p = —5.3, —0.1 — p = —4.1, 


so that convergence depends on the rates at which 


2) daa E 
& 4.1 
get small. Clearly the second one is the more important. 


30.3.5 Summary of Section 30.3 


In this section we defined the terms 


direct iteration with A and with A — pl (page C20) 
inverse iteration with A and with A — pI (page C27) 


We introduced the notation 


m(z) (page C22) 


Techniques 


1. Use direct iteration, with the matrix A — pl, to converge to one 
of the extreme eigensolutions. 


2. Find a p which gives most rapid convergence to one of these 
eigensolutions. 


3. Use inverse iteration with A to converge to the eigensolution 
corresponding to the eigenvalue of smallest magnitude. 


4. Use inverse iteration with A — pI to converge to the eigensolution 
whose eigenvalue is closest to a given p. 


5. Use stable Gauss elimination with interchanges for solving linear 
equations in inverse iteration. 
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30.4 THE METHOD OF SIMILARIT 
TRANSFORMATIONS y 


30.4.0 Introduction 


We turn now to a type of method which is suitable when we want all or 
many of the cigensolutions of A instead of just a few. The basis of all 
such methods is the determination of a non-singular matrix P such that 
the eigensolutions of the matrix B = P~!AP can be found quickly and 
easily. A and B are similar matrices and the transformation from A to B 
is a similarity transformation (of matrices) and therefore preserves eigen- 
values. Moreover, we can easily find the eigenvectors of A from those of 
B, since if y, is an eigenvector of B, then 


By, =2,y,, 
so that 


A(Py,) = P(P~*AP)y, = PBy, = A,(Py,) 
and Py, is an eigenvector of A. 


In using the method we do not try to find the matrix of transition P at 
one step; we perform a succession of similarity transformations which 
simplify the given matrix by stages. 


When the matrix A is symmetric we want our transformations to preserve 
the symmetry. This is achieved by taking P to be an orthogonal matrix, 
so that 


Pol =P". 
Then 
B= PTAP, 


which is also symmetric. In the method we shall describe, each step in the 
sequence of transformations corresponds to a rotation in the plane of 
one pair of coordinate axes. It is possible, with an infinite sequence of 
such transformations (which of course in practice is terminated at some 
suitable point!), to transform A to a diagonal matrix D, whose eigenvalues 
are just its diagonal elements. It is more convenient, however, to use a 
finite number of steps which transform A to a triple-diagonal matrix, with 
non-zero entries only on the main diagonal and on each of the adjacent 
sloping lines. It is fairly easy to determine the eigensolutions of a triple- 
diagonal matrix. In the next sub-section we show how to obtain the triple- 
diagonal matrix, and after that we discuss a technique for finding one or 
more of its eigensolutions. 


30.4.1 Similarity Transformation to Triple-diagonal Form 


We propose to reduce the matrix A to a simpler matrix B by a systematic 
replacement by zeros of the elements of A outside the main diagonal and 
the two adjacent sloping lines. As an example of how the procedure works 
we will look at a 4 x 4 case. Let 


a 2 3 ig 1 0 0 0 

Au| 2 922 423 Ar P = 0 cos 6, —sin@, 0 
x âi}, 423 G33 34 ? ! 0 sind, cos 6, 0 

0o o0 0 1 


ają a24 azs lss 


This orthogonal matrix P, represents a rotation in 4-space through an 
angle @ about the origin in the (2, 3) plane. We shall choose the angle 0, 
to make a pair of elements in the transformed matrix A, = PTAP, equal 
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to 0. Carrying out the matrix multiplication, we find that the (1, 3) and 
(3, 1) entries of A, are 


ai = af) = —5,a,2 + 143, 
where s; = sin 0,, cĉ, = cos 0. 
We can make these entries zero by taking 


aiz i a3 
2 1 A 2-1/2" 
(aia + 473)? (afa + aia)" 


c 


[ae aoe 
Ag+ ay 


ZN 


a2 


The matrix 4;, as you can easily verify, looks like 


a, ayy 0 aig 
A U eee 
rT S a at 


D i 
ajg Sy ON Aga 


all the entries in rows and columns 2 and 3 having been changed by the 
transformation. 


The next step is to make a transformation which reduces another pair of 
entries to zero, while leaving untouched the zeros we have already pro- 
duced in the (1,3) and (3, 1) positions of 4,. This can be accomplished 
by a rotation in the (2, 4) plane, which will change the entries only in the 
second and fourth rows and columns. The new matrix of transition is 


10 0 0 
_|0 ¢ 0 =s 
P=lo 0 1 0 


0 s 0 Cz 
and the new symmetric matrix is 
A; =P] AP. 


You can verify, by carrying out the matrix multiplications, that the (1,4) 
and (4, 1) entries of A, will be zero if we take 


(1) 
Qiz aia 


C= » h= > 
2 TARY aE T (OD? + (ay)? 


The new matrix has the appearance 


a, “7 0 0 


We have now finished with the first row and column, and proceed to con- 


sider the reduction to zero of the entries in the (2, 4) and (4, 2) positions 
of A, This is achieved by the transformation 


A3=P3A2Ps, 
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where Pf is a rotation matrix in the (3, 4) plane, with 
2 

= : a 

2 > 3 Fant, aa 
Naa + RPT MaRS + RPT 


C3 


The entries which change are those in rows and columns 3 and 4, but the 
zeros in the first row and column do not change because the new entries 


in these positions are just linear combinations of the old entries, and 
combinations of zeros produce zeros! 


For our matrix of order 4 we have now fini 


shed, the triple-diagonal 
form looking like 


a b 0 0 
b> a by 0 
c= 0 b; as b, 
0 0 bi a 


a notation which we shall find useful in the next sub-section. 


Example 


12 3 4 12 
3 -12 0 3 
4 0 -12 4 
12 3 4 12 


We will reduce A to triple-diagonal form. At the first Stage a,, = 3, 
a3 = 4, so that c, = 2, s, = $ and we have 


0 
-3 

P,= ; 
0 


coon 
Onemo 
-000 


We calculate A, = PT AP,, and obtain the matrix 


12 5 0 12 
5 -12 0 5 
0 0 -12 0 
12 5 0 12 


A, = 


At the next stage we look at the (1, 4) and (4, 1) entries of A,. a‘) =5, 
a\}) = 12, so that 


Q=% 2=1 
and we take 
I 0 0 
04s 0 -ii 
P=lo 0 1 0 
oyoo 5 


We calculate A, = P7 A;P, to find 


12 13 0 0 
13 12 0 5 
0 0 -12 0 
0 5 0 -12 


At the last stage we reduce the entry 5 in the (2, 4), (4, 2) positions to zero. 
aĝ) = 0, a = 5, so that 


c,=0, s,=1 
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and we take 


100 0 
010 0 
P=|0 00 -1 
001 0 


Performing the appropriate matrix multiplications we calculate 
A, = PI A, P, and obtain 


12 13 0 0 
13:412 5 0 
0 5 -12 0 
0 0 0 -12 


4 = 


which is in triple-diagonal form. Notice that this is rather a special triple- 
diagonal form, two of the entries in the sloping lines next to the diagonal 
being zeros. Clearly A = — 12 is an eigenvalue of A,,'and the others are 
the eigenvalues of the leading 3 x 3-submatrix of A, , which is in “ normal” 
triple-diagonal form. 


For the general case we produce the required zeros in the first row and 
column by successive rotations in the planes (2,3), (2,4),.--,(2,”). 
Zeros in the second row and column are produced by rotations in planes 
(3, 4), (3, 5), ..., (3, n); and so on. The number of successive transforma- 
tions is 

(n-—2)+ (2-3) +--+ 1 =2(n—- 1)(n — 2), 
and we have performed the operation 

PPh -1'0 PEPTAP,P2 «++ Py-1Py = C, 
where C is a triple-diagonal matrix and 

N =n —- 1)(n—-2). 


An eigenvalue of C is an eigenvalue of A, and if y is an eigenvector of C, 
the corresponding eigenvector of A is 


x = (P,P2 ++: Py)y. 


When we have computed an eigenvector y of C we compute the corre- 
sponding x by successively forming Pyy, Py-;(Pwy),---,Pi(P2... Pry, 
rather than by first computing the matrix product (P,P, `+- Py). 


Exercise 


2 4 -— 
Reduce the matrix A= | 4 2 - 
N: eh, Rae 


to a triple-diagonal form. 


Solution 


At the first stage we reduce to zero the entries in the (1, 3) and 
(3, 1) positions. We take 


4 4 3 3 


C EATS s= ETLi -3 


and 
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Then 


25 0 
A,=PTAP,=|5 2 3 
03 - 


which is already in triple-diagonal form. 
30.4.2 Eigenvalues of a Symmetric Triple-diagonal Matrix 


We could use any of our iterative methods to find cigensolutions of a 
triple-diagonal matrix, and the large number of zero entries in this matrix 
considerably reduces the amount of computation involved. But symmetric 
matrices have some special properties which give us a much more attrac- 
tive method, at least for the determination of the eigenvalues. Consider 
the symmetric matrix 


The eigenvalues of these submatrices are approximately 
2, {2.303, — 1.303}, {6.425, — 1.306, — 0.119} 
respectively. 
Inspection of these eigenvalues shows that the eigenvalue of (2] separates 


2 1 
those of |; al: 


— 1.303 < 2 < 2.303; 


and the eigenvalues of (i “3 i] separate those of 


— 1.306 < — 1.303 < —0.119 < 2.303 < 6.425 


This result generalizes to give the first important property of symmetric 
matrices; we state it as a theorem but do not prove it. 


Theorem 1 


For a symmetric matrix, the eigenvalues of the principal submatrix of 
order r separate those of the principal submatrix of order r + 1. 


This result leads to a second theorem, which we again state without proof, 


but which is the basis of a very practical method for computing the 
eigenvalues of a triple-diagonal matrix. 


Theorem 2 


If we take a particular value of p in the matrix A — pl, and compute the 
determinants of its principal submatrices, calling these values Nis fas Sz» 
...,f, for a matrix of order n, then the number of agreements in sign 
between successive members of the sequence fo, fi» - --» Ja, where fo=l, 
is equal to the number of eigenvalues of the matrix A greater than the 


number p. 
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Examples 
(i) With p = 0, the successive principal determinants of the matrix 


are 


fi hh 
1 2 -3 1 


There is one agreement in sign between successive members, between 
Jo and f,, and therefore A has just one positive eigenvalue (since p = 0). 
Gi) With p = —1, we find for the matrix A + J the sequence 


Dhh h ' 
1 3 -1 -2 


Here there are two agreements of sign, which tells us that A has two 
eigenvalues greater than — 1. Since there is only one positive eigenv7lue, 
as we have shown, it follows that there is just one eigenvalue between 
—1 and 0. 


The point of Theorem 2 is that it enables us quite quickly to isolate any 
particular root within a particular interval. To locate it more closely we 
would now proceed to use a process of bisection, taking next p = —0.5 
in our example. If in the resulting sequence there are two agreements in 
sign, we deduce that there is one eigenvalue between —0.5 and O (and 
one greater than 0); and if there is only one agreement, we know that 
there is one eigenvalue in the interval (— 1, — 0.5). 


Once we have located an eigenvalue in an interval (a, b), then after k 
bisections we have confined it to an interval of length 2~* (b — a), and 
convergence to an accurate result is quite rapid. On the other hand we 
have many determinants to evaluate, and this is where the triple-diagonal 
form C shows to real advantage. 


What does the trick, as you might expect, is a recurrence relation! For a 
4 x 4-matrix, using our previous notation, we have 
a-p b 0 0 

b a-p b 0 


C=pl= 0 b, a-p b4 
0 0 bg a-p 
The determinants of the principal submatrices are 
fi =a —-p, 


Sr = (az — p)(a, — p) — b3 = (a2 — p)f, — bf 


where fo = 1. Continuing in this way it can be shown, for a general triple- 
diagonal matrix, that 


Sori = (Gra = DS — Bas fr-i (1) 


a second-order recurrence relation. 


Example 

Our orthogonal similarity transformation applied to the matrix 
1 fifa 32 
v2 -J2 -1 V2 

JA aAA 
2 A2 f2 -3 


A 
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produces, with exact arithmetic, the triple-diagonal form 


1-2/2 0 o 
a -2⁄2 0 J 0 
0 V2 4 ~3 
7 0 -4 -3 


If the eigenvalues are 4,, 22,23 and J, in decreasing order, let us find an 
interval which contains only 2,3. 


We saw in sub-section 30.1.1 that no cigenvalue exceeds M(C) in magni- 
tude. Here M(C) = 5 (from its last row), so that all cigenvalues are in the 


interval [—5, 5]. Let us bisect this and take p=0. Using the recurrence 
relation, we find 


fof hh fs ts 

l 1 -8 2 45 
There are two agreements in sign and therefore two positive cigenvalues, 
which must be 2, and 2,. The required A, is therefore in the interval 


[—5, 0), as also is 2,. Next we bisect this interval, take p = —2.5, and 
compute the sequence 


hhi h fs Ss 
1 3.5 0.75 —5.5 —4.6875 


Three agreements in sign mean that A; is in the interval (-2.5, 0], and 24 
is in the interval [—5, —2.5]. 


Either of these eigenvalues can be located more accurately by successive 
bisections of the relevant intervals. 


Exercise 
The triple-diagonal matrix 


2 -1 0 o0 
-1 2 -1 0 
0 -i 2 -1 
0 0 -1 2 
is positive definite. Find an interval, not exceeding unity in length, which 
contains the smallest eigenvalue. (If any f, =0, its “sign” should be 
taken as opposite to that of f,_,.) 


Solution 


Since C is positive definite all its eigenvalues are positive, by 
Equation (5) of sub-section 30.1.1. Also M(C) = 4; so all eigen- 
values lie in (0, 4). We start with p = 2, and produce the sequence 


to ah h ts Ss 

1 0-1 0 1 
Noting the comment about the “sign” of zero, we find signs 

+--++ 
and the existence of two agreements shows that the two smallest 
roots lie in (0, 2]. 
Bisecting (0, 2), and taking p = 1, we get 

fo hh ty h 

i 1 0 -1 -!I 
Here there are three agreements in sign, and our required smallest 
root lies in the interval (0. 1]. 


30.4.3 Eigenvectors of a Triple-diagonal Matrix 


Once we have found an eigenvalue we can in theory compute the corre- 
sponding eigenvector by solving the linear homogeneous equations 


(C —ANx =0. 
For a matrix of order 3, for example, this system has the form 
(a, — A)x, + b, x2 =0 
by x, + (a, — Ax. + b3x3 =0 (1) 
b,x + (a; — A)x; =0- 
where x = [x, x2 x3]". ; 
Let us first consider the exact solution of a system such as (1) assuming 
that our computed eigenvalue is an exact eigenvalue. We assume that 
ba, b, are both non-zero and look for a solution scaled so that x, = 1. 
Then we can calculate x, from the first equation and x, from the second, 
and the last will then be satisfied automatically since the matrix C — AI 
is singular, so that its rows are linearly dependent. It turns out that x,, x2, 


x, can be expressed in terms of the leading determinants fo , f» /2, for our 
method gives 


x=l=f 

xı = —(a, — a/b, = —filbz 

x, = —(02.%, + (az — A)x2)/b3 (2) 
= ~(bsfo-@- DÉ), 
=fil(b2 by) 


by the recurrence relation (1) of the preceding sub-section. For an n x n- 
triple-diagonal matrix, we obtain by the same method the result 


x, = (1 'F,-1l(b2 bs --- b,) (f=1,...,n) (3) 


Although correct in theory, the method we have just described turns out 
to be quite unsatisfactory for practical computations with inexact arith- 
metic; it suffers from acute induced instability. An example will reveal 
the dangers of this approach. 


Example 


The 21 x 21-triple-diagonal matrix 


10 1 
1 9 1 
1 8 1 
C= i 
1-9 ı 
1 -10 


has an eigenvalue A ~ 10.746 194 183, correct to 1t significant figures. If we 
take the extremely good approximation A = 10.746 194 2, correct to nine 
figures, and denote the corresponding eigenvector by [x, ... X21)", we 
can take x, = 1 and compute the numbers x3, x3,..., X2, in succession 
from the first 20 equations in the system (C — A/)x = 0. The result gives 
a vector some of whose entries are approximately 


X X4 x7 Xn Xis X21 
1 0.09 0.0005 0.036 0.77x 10° 0.19 x 10! 
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The Corresponding entries of the correct cigenvector are completely 
different; they are: 


1 0.09 0.0005 7x1078 2x 107? 0.7.x 10749 


and we have a catastrophic case of induced instability! (This example 


is taken from: J. H. Wilkinson, The Algebraic Eigenvalue Problem 
(Clarendon Press, Oxford 1965).] 


We can explain this result using backward error analysis, by finding the 
equations which our computed “eigenvector” satisfy exactly. Assuming 
no mistakes in the arithmetic, other than a small error in our approximate 
eigenvalue, which we will call 1, we have satisfied exactly the first n — 1 
of the equations (C — A/)x = 0. The last equation is obviously not satis- 
fied, since otherwise both the assumed eigenvalue and the computed 
eigenvector would be exact. Suppose that substitution in the last equation 
produces the number m. Then we have satisfied exactly the equations 


(C — Ax = me, (4) 


where e, is the last column of the unit matrix. What is the solution of 
these equations? Writing 


x=} ax, (5) 


where x,,...,X, are the normalized eigenvectors of C, with eigenvalues 
Ay, «++, Àn, We find that (4) becomes 


È a,(2, — Ax, = me, 
rel 


Taking the inner product with any eigenvector x, and using 


bg giy 
X,’ Xs = XX, = Sess 


we obtain 
(A, — Aa, = me, ' x, (s=1,...,”) 
so that (5) (with s used in place of r) becomes 
_ me," X, 
x= Lao X; (6) 


If Ā is very close to a particular eigenvalue À}, and me, ' xX, is not very 
small, then the term with s = k in (6) will dominate, and our calculated 
vector x is very nearly a multiple of the eigenvector x, which we are trying 
to compute. In our particular example, however, me, - x, is very small, 
since the nth entry in the vector x, is 0.7 x 107'%. Our computed “ eigen- 
vector” is therefore far from being proportional to x,. Such a possibility 
might appear to be unlikely to occur in practice, but experience shows 
that it is in fact quite common with triple-diagonal matrices obtained 
(from symmetric matrices) by orthogonal similarity transformations. 


Now we notice, from Equation (4), that out process is equivalent to one 
step of inverse iteration with a particular starting vector y = me,, and 
it has failed because the matrix C — 2/ is nearly singular and because the 
starting vector is nearly orthogonal to the required eigenvector. Our 
analysis shows that these two causes act together to produce the failure. 
If we remove one of them, by using a starting vector which is not nearly 
orthogonal to the required eigenvector, then the process will succeed and 
usually in one step. This answers a question raised near the end of sub- 
section 30.3.4 about the possible instability of inverse iteration applied 
with an arbitrary starting vector and with the matrix A — pl, when p is 
accidentally close to an eigenvalue. We shall not get instability because 
we shall use several steps of inverse iteration, and rounding errors will 
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guarantee that for some r successive vectors y in the iterative process 
are not orthogonal to the required eigenvector. 


The only danger is with our triple-diagonal C — 1/ in which J is com- 
puted separately and we now seek the corresponding eigenvector. Our 
obvious method is implicitly equivalent to one step only of inverse itera- 
tion, and it is clear that what we have to do, to succeed in just one step, 
is to do this inverse iteration deliberately and properly. In other words we 
solve the equations 


(Cc - ANx = y, 


with Gauss elimination with interchanges, and with a selected y which 
is not nearly orthogonal to the required eigenvector. 


A good choice of y is not easy to guarantee absolutely, but the following 
method has rarely (say once in 10* problems) failed in practice. We 
reduce C — UJ to triangular form U with our stable method of elimination 
with interchanges, so that the equation left for solution is 


Ux = 0. 


Now we are not particularly interested in the precise values of the entries 
of c, but merely that they are obtained from a truly arbitrary y‘. So 
we now take 


cOT (11... 1), 


which corresponds to some particular y, and it would be rather surpris- 
ing if this y were nearly orthogonal to the required eigenvector. This, 
as we have said, almost always succeeds in one step, but for safety (and to 
eliminate the 1 in 10* failure rate), we can quite easily perform just one 
more step of inverse iteration and be virtually certain that the results 
merely confirm the accuracy of the first step! 


Exercise 


In Equation (1) of this sub-section we could work “backwards”, taking 
xX, = 1, computing x, from the last equation and then x, from the second. 
The first equation is then automatically satisfied if A is an exact eigenvalue 
and we perform all subsequent arithmetic exactly. What happens if we 
use this method on the example of this sub-section and with which we 
failed catastrophically with the “ forwards” method? 


In fact, using precisely the same A, working backwards and making no 
mistakes in the arithmetic we find a vector which, suitably scaled, agrees 
with the correct eigenvector to about nine figures! Can you explain this 
fact? ; 


(N.B. This is an accident; the backwards method may fail in other 
examples, just as catastrophically -as the forwards method failed in this 
example!) 


Solution 


In the backwards method all that happens in the previous theory 
is that me, is replaced everywhere by m'e,, where m’ is some num- 
ber and e, is the first column of the unit matrix. Then, in Equation 
(6) of this sub-section 


_ amex, 
ge ae a 


and m'e, ` x, is of reasonable size, since the first entry in x, is 1. 


So our computed x is very nearly a multiple of the Tequired 
eigenvector x, whose eigenvalue 2, is very close to 1. 


30.4.4 Error Analysis of the Similarity Transformation 
Method 


For full analysis of induced instability we should consider: 


(i) the similarity transformation from A to C; 
(ii) computation of the cigenvalues of Cc; 
(iii) the computation of the eigenvectors of C: 
(iv) the “recovery” of the eigenvectors of A. 


The analysis is rather long, so that we shall mention only the main points 
without going into the details. 


In Unit 8 we used a method of backward error analysis to show that the 
computed solution X of Ax = b was the exact solution of 


(A + 6A)K = b + ôb, 


where 5A and ôb have small entries. In the same way, we can show here 
for item (i) that the computed matrix Ĉ, though possibly quite different 
from the exact matrix C computed (if this were possible) with exact arith- 
metic, is nevertheless the exact similarity transformation of some matrix 
A + ôA, where the entries in the matrix 5A are very small in relation to 
those of A. 


There are various contributory factors to this result, but the most import- 
ant is that the successive computed matrices A,, Az,... have entries 
which do not increase significantly in magnitude in relation to those of A. 
In Unit 8 we related the instability of Gauss elimination without inter- 
changes, for solving linear equations, to the growth of the successive 
matrices produced in the reduction of A to upper triangular U. The use 
of interchanges, giving multiples no greater than 1 in magnitude, limited 
this growth to such an extent that we maintained good stability, and the 
computed U was the U which would be obtained by exact arithmetic from 
an original A + 6A with small 5A. 


It is the choice of our particularly simple orthogonal matrices P which 
produces the analogous effect in our transformation of A to the triple- 
diagonal C. 


Item (iv), the recovery of the eigenvectors of A, turns out to be stable for 
quite similar reasons, since we obtain the eigenvectors of A by premulti- 
plying those of C by a product of these simple orthogonal matrices of 
type P. 


For item (ii), the computation of the eigenvalues of C, we can show that 
the use of the recurrence relation for computing 


Siri = ra: — Se — aries 

actually produces an f,,, which satisfies exactly F 
feri = G41 — PS - O fr- 

where al,, and bj+, differ by very small amounts from a,,4, and 6,4; 


respectively. This means that we are finding the exact eigenvalues of a 
matrix differing only slightly from the original triple-diagonal matrix C 


In Items (i), (ii) and (iv), then, there is virtually no induced instability; 
the cigenvalues are always produced with great accuracy, and the eigen- 
vectors of A have little more error than those of C. It is only in item (iii), 
the computation of the eigenvectors of C, that we have to be careful. We 
have already seen in the preceding sub-section how one obvious method 
of calculating the eigenvectors may suffer from induced instability, but 
we have also seen how this instability can be avoided by using inverse 
iteration deliberately and properly. When we use this method of finding 
the eigenvectors of C the entire calculation is therefore free from induced 
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instability, the computed results being the exact results for a slightly differ- 
ent matrix A + 5A with a “small” 5A. Our work on inherent instability 
then guarantees that our computed cigenvalues are very accurate, and the 
cigenvectors will be as good as we can expect in virtue of the possible 
ill-conditioning of eigenvectors corresponding to nearly equal eigenvalues. 


30.4.5 Summary of Section 30.4 


In this section we defined the terms 


similarity transformation (page C31) 
triple-diagonal matrix (page C31) 
leading or principal submatrix (page C35) 
separation of eigenvalues (page C35) 
bisection process (page C36) 


We introduced the notation 


P, (page C31) 
‘A; (page C31) 
5,,¢C, (page C32) 


Cc (page C33) 
tt (page C35) 
Theorems 


1. (page C35) 
For a symmetric matrix, the eigenvalues of the principal submatrix of 
order r separate those of the principal submatrix of order r + 1. 


2. (page C35) 

If the determinants of successive principal submatrices of the matrix 
A — pl are fis f2» ---» Jfa, with fo = 1, then the number of eigenvalues of A 
which are greater than p is equal to the number of agreements in sign in 
successive members of the sequence fo, f,,...,fa- 


Techniques 


1. Find orthogonal plane-rotation matrices to reduce to zero particular 
entries of a matrix by a similarity transformation. 


2. Transform a symmetric matrix to symmetric triple-diagonal form. 


3. Compute principal determinants of a triple-diagonal matrix using a 
recurrence relation, and find the number of agreements in sign in 
successive members of the sequence. 


4. Use this result to determine how many eigenvalues are greater than 
some number, and then to bracket a required eigenvalue in intervals 
of decreasing size. 


5. Avoid the induced instability of an obvious method for computing 
eigenvectors by using one or two steps of inverse iteration instead. 


Formula 


The principal determinants f;, f2, f3, ..., of a triple-diagonal matrix 


a,—p b 0 
b2 a,—p bı 0 
0 b3 a3—p 


satisfy 
Sori = (Gras — PS, — Bai r-1 
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30.5 ERROR AND CORRECTION OF 
APPROXIMATE SOLUTIONS (OPTIONAL) 


30.5.0 Introduction 


The method of backward error analysis, as we know, measures the stab- 
ility of our method but not the accuracy of the computed solution. In 
this section we shall see how to assess the error of an approximate eigen- 


solution, and to improve the eigenvalue, with a relatively small amount 
of extra work. 


In Unit 8 we looked at the similar problem relevant to algebraic equations 
Ax = b, and we observed that the size of the residual vector r = b — AX, 


for an approximate solution X, gave little or no information about the 
error X — x of this approximation. For the symmetric eigenvalue problem, 
however, we shall see that the residual vector turns out to be both useful 
and informative. 


30.5.1 Error of an Approximate Eigensolution 


Suppose that 4 and x are approximations to an eigensolution 2,, x;. 
Assuming x to be normalized, we compute the residual vector 


r= Ax — dx. 


Then, if we define ¢ by 


e= Jen), 
we shall prove that 
Ja, -—Al <e (1) 
To prove this we write, as usual, 
XOX, $5 +4,X, (2) 


where the eigenvectors x,,..., X, are normalized. Then since x, ` x, = 5), 
we have 


l=x:x=a} + +03 (3) 
The corresponding expression for r is 
r = (4 — ADX = Ay — Day +70 + (An — Dan Xn 
from which we obtain 
eè =r: r= (À — Aa t + (An — A) 7a (4) 
> (A, — Aa? + +++ + 22) 


assuming that A is closer to A, than to any other eigenvalue. Using (3) 
and rearranging, we obtain the result (1) which we wanted to prove. 


This is clearly a very useful and easily computed upper bound for the error 
of our computed eigenvalue, and we notice that it does not depend on 
the properties of any other eigensolution. We shall find that the corres- 
ponding result for the eigenvector depends on the nearness of another 
eigenvalue to 2,, and this is hardly surprising when we recall the possibili- 
ties of inherent instability of an eigenvector. 


To estimate the error x — x, in the approximate eigenvector x, we use (2) 
and write 


x — x, = (&; — 1)x, +02. X2 + °° H On Xas 
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so that 
(x — X,)°(x— x) = (a, - 1)? +a? +- +a? 
FAUS 6) 
by Equation (3). 


Our problem is therefore reduced to an estimation of a,. For this purpose 
we define a to be the difference between the approximate eigenvalue 2 
and the nearest eigenvalue other than 2,, so that 


|4, - A] 2a, s=2,3,...,m. 
Then (4) gives 
e = (Ay — APa? + (Az — Aad + e + (A, — A)? 
Baad +--+ +05) n 
=a*(1 — a3). (6) 


So 
2 


E 
asl- 


and if we neglect fourth and higher powers of e and assume, reasonably 
enough, that @, is positive, this gives simply 


g? 
2(1 -«,) < w 
and (5) finally gives the required result 
(x) a-si G 
This says that the sum of the squares of the entries in x — x, is at most 


2 
8 PERAN i ; A 2 
w and the individual entries, the errors in those of x in relation to those 


e. ; 
of x,, can hardly exceed z in magnitude. 


Example 
In the Example of sub-section 30.3.2 we found the approximations 
y = [0.731 0.233 1.000)", 4 = 6.43 


for the largest eigensolution of the matrix 


2 -3 
A=ļi -iI 1}. 
3 1 4 


Let us find upper bounds: for the errors of these results. 


The required normalized approximation x is y/k, where k = (y- y)'/?, 
but the arithmetic is simplified (and more accurate!) if we defer the norm- 
alization. 


Computation gives 


k? =y-y = 1.588 650, 
and 


r = (4 —Alx = (A — ANy/k 


=; [—0.005 33 —0.00019 —0.00400]" 


so that 
2 _ „ . p _ (0.005 33)? + (0.000 19)? + (0.004 00)? 
e =r'r= E 
_ 0.000 044 445 
1.588650 
= 0.000 027 98 
to four significant figures, giving 
e = 0.0053. 


By Equation (1) the error in the computed eigenvalue is then less than 
0.0053 in magnitude. Moreover, since the eigenvalues of A are approxi- 
mately 6.4, —1.3 and —0.1, the next nearest eigenvalue is at a distance 
of about 6.5 from A,, So that a ~ 6.5 and the error in x, the normalized 


approximate eigenvector, cannot (by (7)) exceed e/a ~ 0.0008 in any 
entry. 


Accurate computation shows that the error in the eigenvalue is 0.0050, 
and the maximum error in the eigenvector is 0.0003 in its first entry. Our 
error bounds, 0.0053 and 0.0008, do not therefore greatly overestimate 
the true errors. 

Exercise 

Repeat the process of the example with the approximation 


y=(0.7 0.2 1.07, A=6. 


Solution 
We find . 
(A — 6Dy = [0.4 0.3 0.3)" 
and 
yy = 1.53, 
so that 


0.34 
2m = 0,222 
e =s 


and 
e = 0.47, approximately. 


So the error in the eigenvalue is at most 0.47, and in the nor- 
malized eigenvector at most 0.08 in any component. 


30.5.2 Correction of an Approximate Eigensolution 


Finally, we examine the possibility of improving our approximation with 
little extra work. This we can do very easily for the eigenvalue, the more 
accurate estimate corresponding to our approximate eigenvector y being 
the Rayleigh quotient 


J: 
R= y = = xTAx 
where 
y 
=y y” 
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is the corresponding normalized approximate eigenvector. The reason 
why the Rayleigh quotient is a good estimate can be seen by writing 


X= QjX, HAX, +°°°+4,X, 
so that 

Ag = apd, 03d, ++ + Rd, (1) 
by Formula (5) of sub-section 30.1.1. It follows, since 

aft tay = 1, 
that 

Ag — Ai = (at Ar Hee + andy) — (a1 + 0+ + aA 

= a3(Ay — Ay) +++ + F(A, — Ay) 


It follows that if « denotes the error in x as an approximation to x,, 
measured by the largest of the numbers |a.|,....,|a,], then the error 
Ag — A, in the improved approximate eigenvalue A, is roughly a’, and is 
likely to be much smaller than the error of the approximate eigenvector 
it was obtained from. 


We can even find a bound for the error [Ag — ,| of the Rayleigh estimate. 
Using (1) and the normalization condition 


a+ +o =1, 
we can write 
Aga? + °° + a?) = Ayo? +-++ + Aa? 
so that 
Qr — Ay)ot = (42 — Aaah + +++ + (A, — Ande? 
Q2 — Ax)*a3 A, — r)a 
we a a cg a 
à2— År 4, — Àr 
It follows that 
Qaia... On = Aan 
lå2 — årl [an — Arl 


1 
<5 (a — Anya +--+ + Ay — Aado] 


làr — Ay lat < 


where a is the smallest of the denominators. Defining 
r= (A —A,)x 
= Ar — Ander X, +277 + (An — Aga Xn > 


we obtain 


1 1 
lår — Ayla? <r =e? 
a a 


so that 


e? 
lår -ål < m 
5 2 
<——> V f sub-secti 
<a AA y (6) of sub-section 30.5.1. 
A reasonable knowledge of a thus gives a very useful estimate, approxi- 
mately ¢?/a, for the error of the Rayleigh quotient. 


Unfortunately we have no simple method for correcting or improving 
our approximate vector. We were able to improve the eigenvalue by 
evaluating the Rayleigh estimate because the error in our approximation 
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did not depend on any other eigensolution; but we have no simple methods 
for correcting the Rayleigh estimate, or of obtaining any simple correction 
to the approximate eigenvector, because the error bounds for these approx- 
imations do need a knowledge of the other eigensolutions. Knowing 
approximations to all the eigensolutions we can in fact correct the approx- 


imate eigenvectors but this is a much more elaborate operation which 
we shall not consider here. 


Example 


Let us find (i) the Rayleigh estimate for the approximate eigenvector 
y = (0.731 0.233 1.000]7 
for the matrix 


used in the Example of sub-section 30.5.1, and (ii) a bound for its error. 
We find 
Ay = [4.695 1.498 6.426)" 
and with 
yTy = 1.588 650 
we get 


y7Ay 10.207 079 ; 
i = 1.588 650 7 7 tely. 
An= Sry = 7598 650 = 425 001 7 approximately. 


For the residual vector we find 


r= (y+ y) (Ay —Apy) ; 
= (y+ y)~/?[—0.001 68 0.00097 0.001 00]7 


(to five decimals), so that 
è ge 
€? = 0.000 003 00, — ~ — ~ 0.000 000 5. 
a 6.5 


This is an error bound for the Rayleigh quotient, its actual error being 
about 0.000 000 4. We note the great improvement, in our estimate of the 
eigenvalue, which we obtain by computing the Rayleigh quotient. 


Exercises 


1. If y is an approximate eigenvector, À an approximate eigenvalue, and 


r = Ay — 2y, show that the Rayleigh quotient gives a correction to 2 
T, 


r 
of amount AR 
Af 3 
2. Using the method of Exercise 1, and the result of the exercise of the 


previous sub-section, show that, for the vector 
y=[0.7 0.2 1.0)’, 
the Rayleigh quotient gives Ag = 6.42, approximately. 


Solutions 


1. Since 


Ay — Ay =r, 
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and 


y"Ay 
A= 
oyy 
we have 
y@yt+n aI 
'R T, 


2. The exercise of the previous sub-section gives 
r=(04 03 0.3)”, y7y = 1.53, and 2 = 6, 
so that 


(0.4)(0.7) + (0.3X0.2) + (0.3X1.0) 
mót ( 1.53 ) SONN 


to two decimal places. This has a maximuni error of little 
more than 0.005, although the vector is correct to only one 
significant figure. 


30.5.3 Summary of Section 30.5 


In this section we defined the term 
Rayleigh quotient (page C45) 
and introduced the notation 
An (page C45) 


Techniques 
1. Calculate a bound for the error of an approximate eigenvalue by 
computing the residual vector of an approximate eigensolution. 


2. With a knowledge of the next nearest eigenvalue, calculate a bound 
for the error of the approximate eigenvector. 


3. Improve the approximate eigenvalue by calculating the Rayleigh 
quotient, and determine an upper bound for the error of this approxi- 
mation to the eigenvalue. 


LM 30.6 
30.6 SUMMARY OF THE UNIT 


In the first section we established some properties of eigenvalues and 
eigenvectors, for symmetric matrices, needed in the subsequent numerical 
work. These included measures for the “magnitudes” of matrices and 
vectors. We also observed, and illustrated, the lack of economy and stab- 
ility in the use of the characteristic equation for computing eigenvalues. 


The second section discussed the problem of inherent instability, and we 
showed that the eigenvalues are always “ well-conditioned ” in relation to 
small perturbations in the entries of the matrix, and that the eigenvectors 


are ill-conditioned only when associated with eigenvalues which are nearly 
equal. 


In the third section we discussed a method of direct iteration with the 
matrix A — pl, choosing p to give convergence to one or other of the 
extreme eigensolutions, the only possibilities with this method. Then we 
showed that inverse iteration with A — pI (direct iteration with (A — pI)7') 
will converge to the eigensolution whose eigenvalue is nearest to P, so that 
by suitable choice of p we can find any eigensolution. The iterative methods 
were expressed in practical computational form, and we illustrated the 
method of Gauss elimination with interchanges for solving the linear 
equations involved in inverse iteration. We also showed that all these 
methods are virtually free from induced instability. 


In the fourth section we showed how we could transform the given matrix 
in a finite sequence of orthogonal similarity transformations into a sim- 
pler form, the triple-diagonal form which is very suitable for computation. 
We showed how to compute determinants of successive principal sub- 
matrices of C — pI by using a simple recurrence relation, and observed 
that by inspecting the signs of members of this sequence we could find 
the number of eigenvalues exceeding the number p. A systematic bisection 
process then located an eigenvalue quickly and accurately. The obvious 
method of computation of the eigenvector, however, we found to exhibit 
significant induced instability, and replaced this by inverse iteration which 
always converged in one or at most two steps. With this refinement all 
our methods were free from induced instability. 


In the final (optional) section we found bounds for the error in an approxi- 
mate eigensolution, and showed how the computation of the Rayleigh 
quotient, an economic operation, gave a much better estimate for the 
eigenvalue. 


Definitions 

eigensystem (page C5) oom 
modal matrix (page C7) = 
scaling of matrices (page C9) tx 
perturbation of a matrix (page C13) E% 
first-order perturbation of an eigenvalue (page C13) i 
first-order perturbation of an eigenvector (page C13) ats 
direct iteration with A and with A — pI (page C20) Ee 
inverse iteration with A and with A — pI (page C27) es 
similarity transformation (page C31) 7 
triple-diagonal matrix (page C31) * 
leading or principal submatrix (page C35) i 
separation of eigenvalues (page C35) g 
bisection process (page C36) 

Rayleigh quotient (page C45) 
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Theorems 


1. (page C35) 
For a symmetric matrix, the eigenvalues of the principal submatrix of 
order r separate those of the principal submatrix of order r + 1. 


2. (page C35) 

If the determinants of successive principal submatrices of the matrix 
A — pl are fi, fz, -- -» fa, With fo = 1, then the number of eigenvalues of A 
which are greater than p is equal to the number of agreements in sign in 
successive members of the sequence fo, fi, -- -s fa- 


Techniques 
1. Scale a matrix so that the magnitude of the largest entry of the new 
matrix is 1. 


2. Use direct iteration, with the matrix A — p1, to converge to one of 
the extreme eigensolutions. ‘ 


3. Find a p which gives most rapid convergence to one of these eigen- 
solutions. 


4. Use inverse iteration with A to converge to the eigensolution corre- 
sponding to the eigenvalue of smallest magnitude. 


5. Use inverse iteration with A — pI to converge to the eigensolution 
whose eigenvalue is closest to a given p. 


6. Use stable Gauss elimination with interchanges for solving linear 
equations in inverse iteration. 


7. Find orthogonal plane-rotation matrices to reduce to zero particular 
entries of a matrix by a similarity transformation. 


Transform a symmetric matrix to symmetric triple-diagonal form. 


Compute principal determinants of a triple-diagonal matrix using a 
recurrence relation, and find the number of agreements in sign in 
successive members of the sequence. 


10. Use this result to determine how many eigenvalues are greater than 
some number, and then to bracket a required eigenvalue in intervals 
of decreasing size. 


11. Avoid the induced instability of an obvious method for computing 
eigenvectors by using one or two steps of inverse iteration instead. 


12. Calculate a bound for the error of an approximate eigenvalue by 
computing the residual vector of an approximate eigensolution. 


13. With a knowledge of the next nearest eigenvalue, calculate a bound 
for the error of the approximate eigenvector. 


14. Improve the approximate eigenvalue by calculating the Rayleigh 
quotient, and determine an upper bound for the error of this approx- 
imation to the eigenvalue. 


Notation 
Amaz (page C8) 
M(x) (page C8) 
M(A) (page C8) 
A(e) (page C13) 
x,(e) (page C13) 
m(z) (page C22) 


P, (page C31) 
A, (page C31) 
Sr, (page C32) 
C (page C33) 
I (page C35) 
AR (page C45) 
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30.7 SELF-ASSESSMENT 


Self-assessment Test 


This Self-assessment Test is designed to help you test your understanding 
of the unit. It can also be used, together with the summary of the unit, 
for revision. The answers to these questions will be found on the next 


non-facing page. We suggest that you complete the whole test before 
looking at the answers. 


(In these questions all the matrices are symmetric.) 


1. (i) For any vector y we know that 
YAY < YY jhn 
where 2,,., is the eigenvalue of largest magnitude of a matrix A. 
For what vector y do we have equality in this expression? 
(ii) (optional) 
If y is near to some particular eigenvector of A, what can we 


YAY, 
yy 


say about the number 


(iii) Define the number M(A) for a matrix A, and say what value 
M(A) has for the matrix 


1 -2 
a[i -A 
Verify that 
M(A) > |A;| (= 1, 2) 


where A, and A, are the eigenvalues of A. 


2. (i) A matrix A, has eigenvalues 3.1, —0.6 and —0.8, and a matrix 
A, with the same eigenvectors as A,, has eigenvalues 2.9, —3.8, 
and —4.8. For which of these matrices are we more likely to 
have significant ill-conditioning in (a) the eigenvalues, and 
(b) the eigenvectors, in respect of the same small perturbations 
in the entries of the matrices? 


(ii) Which of the eigenvectors are likely to exhibit the greater degree 
of ill-conditioning in 2(i), and what is the approximate ratio of 
the magnitudes of this ill-conditioning for the matrices A, and 
Az? 


3. We propose to compute the eigensolutions of the matrices A, and A, 
in Question 2(i) by an iterative method. 


(i) To what eigenvalues shall we converge using direct iteration 
with A, — pl and A, — pl for 


(a) p=0, 
(bl) pesi 
(c) p=50? 
(ii) What values of p would we use to obtain most rapid convergence 
to 


(a) the eigenvalue 3.1 of A,, 
(b) the eigenvalue —4.8 of A3, 
(c) the eigenvalue —3.8 of A,? 


(üi) If we use inverse iteration with A, — pI and A, — pl, to what 
eigenvalues shall we converge with p=1, and how can we 
estimate the rate of convergence? 
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4. (i) Reduce to triple-diagonal form C, by means of an orthogonal 
similarity transformation, the matrix 


1 2 2 
A=|/2 -/2 -1 
v2 -1 V2 


Gi) An eigenvector of C is 
x= [- Z 0 ni 
Write down the corresponding normalized eigenvector of A. 
(iii) Find the number of positive eigenvalues of A. 
5. The eigenvalues of the triple-diagonal matrix 
1 2 0 
C=|2 -1 2 
0/2 1 
are arranged in order 2, > 4, > A;. 


(i) All the eigenvalues lie in the interval (—a, a). What is the 
smallest value of a that you can obtain almost by inspection 
from the matrix C? 


(ü) By taking a = 6, for simplicity, find an interval which contains 
only the eigenvalue 2, . 
6. (Optional) 
For the matrix of Question 5 we have computed the approximate 
eigensolution (with normalized eigenvector) 


(i) Show that the error of this approximate eigenvalue does not 
exceed (9 — 6,/2)!/2. (This is just larger than 0.7, and the cor- 
Tect eigenvalue is Ji ~ 2.65.) 

Gi) Find the Rayleigh quotient for this vector, and observe its great 
accuracy as an approximation to the eigenvalue (the approxi- 


mate eigenvector has only the first figure accurate in all three 
components). 
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Solutions to Self-assessment Test 


1, 


3. 


G) 
Gi) 


Gii) 


Gi) 


@ 


Gi) 


ii) 


G) 


The required vector is any multiple of the eigenvector corres- 
ponding to the eigenvalue Amax- 

The number y7Ay/y’y is the Rayleigh quotient. If y is close to 
some eigenvector then the Rayleigh quotient is even closer to 
the corresponding eigenvalue. 

M(A) is the maximum row sum of absolute values of the entries 
of A. For the given matrix M(A) = 6. The eigenvalues satisfy 


(1 — A)(4. — 4) — 2? =0, or 4? — 5A =0. 
Then 4, = 0, 2, = 5, and 
M(A) > [àl (= 1,2) 


(a) Neither matrix exhibits significant ill-conditioning in the 
eigenvalues, and the “first-order” perturbations in the 
eigenvalues corresponding to the same eigenvector have 
the same value for each matrix. 

(b) Two of the eigenvectors of A, are more likely to be ill- 
conditioned than any eigenvectors of A, in virtue of the 
near equality of the eigenvalues —0.6 and —0.8 of A). 
The eigenvalues of A, are well separated and its eigen- 
vectors therefore well-conditioned. 

Each of the eigenvectors corresponding to the eigenvalues —0.6 

and —0.8 of A, will exhibit ill-conditioning, by a factor pro- 

portional to 1/(—0.6 + 0.8) = 5. The corresponding number 
for the matrix A, is 1, so the relevant eigenvectors of A, are 
likely to be “5 times as ill-conditioned” as the corresponding 

eigenvectors of A. 


(a) With p =0, we converge with direct iteration with A, to 
the eigenvalue 3.1, and with A, to the eigenvalue —4.8, 
those farthest from the origin. 

(b) With p=-—1, we converge with direct iteration with 
A, — pI to the eigenvalue 3.1 and with A, — pI to the 
eigenvalue 2.9, those for which |à, — p| is largest. 

(c) With p = 50, we get convergence, for the same reason, to 
the eigenvalues —0.8 and —4.8. 

(a) and (b) The required values of p are 


4(—0.6 — 0.8) = —0.7, and 4(2.9 — 3.8) = —0.45. 


(c) No value of p will give convergence to the eigenvalue — 3.8, 
the “non-extreme” eigenvalue of A3. 

With p = 1 we converge to the eigenvalue for which |4, — p| 

is smallest, that is the eigenvalues —0.6 of A, and 2.9 of A3. 

The rates of convergence depend on the ratios of 


[Gx — PYG, — p), 


where A, is nearest to p, and 4, is another eigenvalue. For the 
first result of (iii) we converge to the eigenvalue —0.6 at a rate 


., (1.6\% 1.6\* 
at which (3) and (5) get small as s increases and for the 


s s 
second result the corresponding numbers are (3) and (3) : 


The required orthogonal matrix is a plane-rotation in the (2, 3) 
plane through angle 0, where 


cos 8 = a3, /(aj, + 43,)*”?, 


sin 6 = ay,/(a3, + a3,)'?, 
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(ii) 


(iii) 


G) 


Gi) 


giving cos 0 = sin 0 = A 
2 


1 0 0 


1 
Then P= g V2 
1 


vale fi 


and by computation we find that the required triple-diagonal 
form (only one rotation being needed) is 


1 2 0 
PTAP=C=|2 -1 ./2 
0/2 #1 


The eigenvector of A corresponding to an eigenvector x, of C 
is Px,, which here is 


1 1 

sas at: a: 
eT o Pt 
v2 V2 | v2 
Oe ell 4 ae 
v2 2 v2 


This eigenvector x, normalized so that x?x = 1, is therefore 


=U 1 I). 


WE 


The number of positive eigenvalues of A is the same as the 
corresponding number for C. Using the relevant recurrence 
relation (or by direct calculation) we find the sequence 


hhh h 
1 1 -5 -7 


for successive determinants f,, f2 and /3 (with fọ = 1) of principal 
submatrices of C — OI = C. There are two agreements in sign 
and therefore two positive eigenvalues (greater than p = 0). 


No eigenvalue exceeds M(C) in magnitude, so that a = M(C) 
is the smallest easily calculable value of a. We find M(C) = 
3+ J2 (from the second row). 


Taking p = 0 (bisecting the interval (— 6, 6)) we find for succes- 
sive principal submatrices of C the sequence 


hhh h 
1 1 -5 -7 


with two agreements in sign. Then there are two eigenvalues 
in the interval (0, 6), which must be A, and 22. We therefore 
take p = 3, and find for C — 37 the sequence 


hhh fs 
1-2 4 -4 


There are no agreements in sign, so that both positive roots are 
in the interval (0, 3]. 
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Gi) 


Now take p = 1.5, and for C — 1.5] we find 
h h h fs 
1 -05 —2.75 2.375 
There is now one agreement in sign, so that A, is in the interval 
(1.5, 3], and A, in the interval (0, 1.5). 


Since the approximate eigenvalue is normalized, the error 
bound for the approximate eigenvalue À is (r-r)'?, where 
r=(C—ADx. We find 


Pe (l-$/2 -2+3/2 ~4+4V/3, 


and a little computation gives r-r = 9 — 6/2. The error bound 


for Xis therefore (9 — 6/2)". 
Since the approximate eigenvector is normalized, the Rayleigh 
quotient is just x7Cx, which is equal to 4 + x7r. We find 


xr = 1/2 — 1) = 0.62, 


so that the Rayleigh quotient is approximately 2.62, agreeing 
to two figures with the correct eigenvalue. 
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